## Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan

library(R2jags)

odd <- function(x) x%%2==1

J <- 4 ## number of issue areas asked about
K <- 5 ## number of endorsement groups (including control)
N <- 1500 ## number of survey respondents
n.province <- 4
n.division <- 16
province <- c(rep(1,n.province), rep(2,n.province), rep(3,n.province), rep(4,n.province))
n.per.division <- rep(round(N/n.division),n.division); n.per.division[n.division] <- N - sum(n.per.division[1:(n.division-1)])

delta.division <- c(-1.00,-0.75,-0.50,-0.25, #province 1
		-0.40,-0.10,0.10,0.40,  #province 2
		-0.75,-0.25,0.25,0.75,  #province 3
		0.00,0.50,0.50,1.00)  #province 4

n.age <- 4
delta.age <- c(0.5, 0.5, 0, -1)
age <- sample(1:4, size=N, replace=T)

delta.female <- 0.5
female <- rbinom(N, 1, p=0.5)

lambda.division <- matrix(c(rep(0,n.division), rep(0, n.division), rep(-delta.division, K-3), delta.division/3), n.division, K)
theta.age <- matrix(c(rep(0,n.age), 0.5,0.25,0,-0.75, 0.25,0,0,-0.25, rep(c(-0.75,-0.25,0.25,0.75),K-3)), n.age, K)
theta.female <- c(0, -0.25, 0.5, rep(0, K-3)) ## length K
omega <- c(0, rep(0.3, K-1))


division <- NULL
for (i in 1:n.division){  division <- c(division, rep(i,n.per.division[i])) }
## Ideal Points

x <- rep(NA, N)
for (i in 1:N){ x[i] <- rnorm(1, mean=(delta.age[age[i]] + delta.female*female[i] + delta.division[division[i]]), sd=1) }
x <- x / sd(x)

beta <- runif(J, 0.5, 2) ## discrimination parameter
alpha <- matrix(0.1,J,4)
for (j in 1:J){  ## variation in thresholds, but distance ensured
	while(alpha[j,2]-alpha[j,1] < 0.5 | alpha[j,3]-alpha[j,2] < 0.5 | alpha[j,4]-alpha[j,3] < 0.5){
		tmp <- rnorm(4, mean=0, sd=1)
		alpha[j,] <- tmp[order(tmp)]
	}}

s <- array(NA, c(N,J,K)) ##Individual level impact of group endorsement k
for (i in 1:N){
	for (j in 1:J){
		for (k in 1:K){
			s[i,j,k] <- rnorm(1, mean=(theta.age[age[i],k] + theta.female[k]*female[i] + lambda.division[division[i],k]), sd=omega[k])
		}}}


##Ordered Responses (5 choices)
Y <- array(NA, c(N,J,K))
for (i in 1:N){
	treat <- ifelse(odd(i),TRUE,FALSE) ## Half of Jake's respondents receive the control, remaining half split among treatments
	for (j in 1:J){
		for (k in 1:K){
			p1 <- pnorm(alpha[j,1] - beta[j]*(x[i] + s[i,j,k]))
			p2 <- pnorm(alpha[j,2] - beta[j]*(x[i] + s[i,j,k])) - (p1)
			p3 <- pnorm(alpha[j,3] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2)
			p4 <- pnorm(alpha[j,4] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2 + p3)
			p5 <- 1 - (p1 + p2 + p3 + p4)
			
			p <- c(p1,p2,p3,p4,p5)
			Y[i,j,k] <- sample(c(1:5), size=1, prob=p) ## Full array of potential outcomes
		}}
	if (treat==FALSE) { Y[i,,2:K] <- NA}
	if (treat==TRUE) {
		Y[i,,1] <- NA 
		for (j in 1:J){ 
			answered <- sample(2:K,1); Y[i,j,-answered] <- NA  
		}
	}}

Y2 <- matrix(NA,N,J)
endorser <- matrix(NA, N,J)
for (i in 1:N){
	for (j in 1:J){
		endorser[i,j] <- (1:5)[!is.na(Y[i,j,])]
		Y2[i,j] <- Y[i,j,endorser[i,j]]
	}}
Y <- Y2


lambda.division.real <- lambda.division
theta.age.real <- theta.age
theta.female.real <- theta.female

#######################################
##  Prep and send data through JAGS  ##
#######################################
data <- list ("N", "J", "K", "Y", "endorser",
		"division", "n.division", 
		"n.age", "age", "female",
		"province", "n.province",
		"lambda.division.real", "theta.age.real", "theta.female.real") 
inits <- function() {list(alpha0=matrix(seq(-2,2,length.out=(4*J)),4,J), beta=runif(J, 0.2, 2), 
			x=rnorm(N), s=matrix(rep(0, N*J), N,J), 
			delta.division=rep(0,n.division), delta.province=rep(0,n.province),
			delta.age.raw=rnorm(n.age), delta.female=rnorm(1),
			xi.age=runif(1),
			theta.age.raw=matrix(c(rep(0,n.age), rnorm(n.age*(K-1), n.age, K)), n.age, K), theta.female=c(0,rnorm(K-1)),
			lambda.division=matrix(c(rep(0,n.division), rep(0,n.division*(K-1))), n.division, K), theta.province=matrix(c(rep(0,n.province), rep(0,n.province*(K-1))), n.province, K),
			xi2.age=c(0,runif(4)),
			omega=c(0,runif(K-1)) )} 
params <- c("lambda.division.error", "theta.age.error", "theta.female.error") 
model <- "ModelSimulationsDivisionWithICovariates.txt" 
fit <- jags(data, inits, params, n.iter=20000, model.file=model)
fit <- autojags(fit, n.iter=20000, n.thin=20, Rhat=1.02, n.update=5)
save.image(file = Outfile)

